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Abstract 

A Feynman-alpha formula has been derived in a two region domain 
pertaining the stochastic differential self-interrogation (DDSI) method 
and the differential die-away method (DDAA). Monte Carlo simula- 
tions have been used to assess the applicability of the variance to 
mean through determination of the physical reaction intensities of the 
physical processes in the two domains. More specifically, the branch- 
ing processes of the neutrons in the two regions are described by the 
Chapman - Kolmogorov equation, including all reaction intensities 
for the various processes, that is used to derive a variance to mean 
relation for the process. The applicability of the Feynman-alpha or 
variance to mean formulae are assessed in DDSI and DDAA of spent 
fuel configurations. 

1 Introduction 

Special nuclear materials have been assessed using the differential die-away 
analysis (DDAA) method [U [2j [3j E3 [5] for several decades. In this determin- 
istic method an initial pulse of fast neutrons is injected to the sample and 
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the time dependence of the detection rate of fast neutrons is used to deter- 
mine the presence of fissile materials such as 235 JJ and 239 Pu. Its application 
is suitable when the fissile material is embedded in a moderating surround- 
ings such that the source neutrons induce thermal fission after having slowed 
down. 

There exists, however, one significant drawback of pulsed measurements, 
that is the necessity of a neutron generator. The stochastic generalization of 
the DDAA method was recently suggested as an alternative method which 
eliminates the need for such an external source, namely the so-called Differ- 
ential Die-away Self-Interrogation(DDSI) technique [B]. The DDSI method 
utilizes the inherent spontaneous neutron emission of the sample. In the ab- 
sence of a trigger signal, the temporal decay of the correlations as a function 
of the time delay between two detections of fast neutrons is used in the DDSI 
method. This corresponds to a Rossi-alpha measurement with two energy 
groups. 

An empirical two-group Rossi-alpha formula, representing the functional 
form of the DDSI formula, was derived in Ref. [6j. This result was subse- 
quently re-derived from first principles with the use of the backward-type 
Kolmogorov equations [7J |8]. Note that, the two-group theory of the Rossi- 
and Feynman alpha formulae is interesting also in other areas than nuclear 
safeguards, such as pulsed and stationary source driven experiments mea- 
suring the reactivity in fast cores of accelerator driven sub-critical systems. 
In such experiments [HI HDl E] it was found that two exponentials appear, 
indicating that the temporal behavior of the fast and thermal neutrons is 
separated, especially in fast reflected cores. This amplified the need for the 
two group versions of the Feynman and Rossi-alpha formulae [8]. 

Energy dependent aspects of neutron counting have been studied pre- 
viously; however, most investigations are limited to formulae valid for one 
infinite region. Note that there exists exceptions, in Ref. [12] a two region 
theory representing the core and a reflector region is discussed. In certain 
situations one cannot disregard the fact that the neutron processes are taking 
place in a non-homogeneous medium, which consists of different regions with 
different characteristics. For instance in the DDAA and DDSI methods, the 
fissile material is surrounded with a moderator. Hence the fission takes place 
in the inner region, whereas the slowing down in the outer one, and there is 
a current of neutrons between the two regions in both directions. 

If the domain is divided into several regions, two or more distinct expo- 
nentials will appear in the Feynman- and Rossi alpha formulae. In this con- 
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tribution, a stochastic theory for a branching process in a one group neutron 
population applied to two domains is studied, based on the previous results 
of Refs El H3] . In particular, we consider a counterpart of the Differential 
Self-interrogation Variance to Mean (DSVM) formula, that is derived by us- 
ing the master equation or Kolmogorov forward approach [T3]. The model 
includes a spontaneous fission source of fast neutrons, absorption in both 
domains, thermal fission, and detection of fast neutrons. Furthermore, we 
have included the intensity of neutron passage from one region to the other. 
However, the present variance to mean formula is energy independent and in 
contrast to the energy-independent one point theory of multiplicity, where 
the only appearing parameter is the first collision probability, we find two 
distinct exponentials appearing. 

The usefulness of the theoretically predicted variance to mean, that is 
based on the observation of the two exponentials depends on the quantitative 
values of the various within- and inter-region neutron reaction intensities, 
which appear as coefficients in the equations. To assess the applicability 
and expected performance of the DSVM method in practical situations, the 
above mentioned reaction intensities were determined from numerical Monte 
Carlo (MCNPX) simulations. Significance of different values of the reaction 
intensities of thermal and fast neutrons in the performance of the method is 
discussed. 

2 A two-point variance to mean formula 

In Ref. JT3], the authors discussed the extension of the variance to mean 
or Feynman-alpha one-group to two-group theory in a stochastic setting de- 
scribed by the Kolmogorov forward approach. Here, we will instead consider 
a one-group theory extended to two domains, also known as two-point the- 
ory. In order to keep the possibility of analytical solutions, we will choose a 
description which still disregards the handling of the spatial transport. The 
model being considered is that of two adjacent homogeneous half-spaces with 
space-independent reaction intensities in both regions. The exchange of the 
neutrons between the two regions will be described by the two passage inten- 
sities into the two different directions. The philosophy behind this model is 
the same as in the two- or many-point theories of coupled reactor cores, first 
introduced by Avery [TJ] where a two-region system is described by point 
kinetics in each region, whereas the interaction between the two regions is 
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described by coupling constants. Obviously, the assumption of a generic pas- 
sage intensities which are the same for all neutrons in the respective regions 
is only an approximation of the real case, and would only correspond to the 
physical situation if both regions were infinite and took up the same space 
by intertwining each other. Nevertheless, the approximation in neglecting 
the space dependence of the transition intensity is not significantly coarser 
than using one single infinite homogeneous model for the evaluation of mea- 
surements performed in finite systems. 

In the model, we have included a Poisson source of fast neutrons described 
by the source strength Si that releases n particles with probability p q (n) at 
an emission event (i.e. spontaneous fission) in region I. Furthermore, we 
assume that the source is switched on at time t — to < 0, although the 
dependence on to will not be denoted. The detection rate of particles is 
included and is denoted by the intensity Xd- We will start by deriving the 
explicit analytical formulae describing all the possible processes in the system. 
The Kolmogorov forward equation will give us a differential equation for the 
probability P(N\, N 2 , Zi,t) for having N\ neutrons in region I, N 2 neutrons 
in region II at time t in the system and having detected Z\ neutrons in 
the interval (0,t). In deriving this differential equation we have summed all 
mutually exclusive events during an infinitesimally small time interval dt and 
we find for the probability P(N 1 ,N 2: Z 1 ,t) the differential equation 
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Here, f 3 k is the probability of having exactly k neutrons produced in an 
induced fission event in region j and Ai and A2 are the decay constants 
(total reaction intensities) for particles in region I and II, respectively. The 
intensities Ai a , A 2a are the absorption (actually, capture) intensities of region 
I and II, while fission resulting from the thermal particles happens with the 
intensity of Aij and A2/. The intensity of the detection of a fast neutron is 
denoted by A^. The total intensities are given by 

Ai = A ia + Ai/ + An + A d , (2) 

and 

A2 = A2a + A2/ + At2- (3) 

Where, An describes the intensity of particles leaving for region II and An 
is the intensity of particles transferring to region I. We derive the equations 
for the factorial moments by using the generating function of the form 

c(x,Y,z,t) = j2Y,J2 xNlYN2ZZip ( N ^ N ^ z ^> ( 4 ) 

Ni N 2 Zi 

and describe the time evolution of the process by a partial differential equa- 
tion in the variables (X, Y, Z) in terms of the generating function as, 

BC BC 
— = (\ la + q \X)\ lf + \ d Z~\ 1 X-\ T1 Y) — 

BC 

+ (X 2a + \2 f q 2 (Y)-\ 2 Y-\ T2 X) — 

+ S l (r(X)-l)G, (5) 

where 

Q ] {X) = zZf k X\ and (6) 

k 

r(X) = ^p,(n)I«. (7) 

n 

Here the superscript j denotes the possibilities of having different material 
compositions in region I and II, respectively. In the numerical analysis we 
will assume that q 2 (Y) vanishes signifying that region I is surrounded by 
moderating material without fissile components. We have used the definition 
of the derivatives of the expressions (p| and (17) as v\ = dq> / dX\x=i and r\ = 
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dr / dX\x=\, which stand for the expectations of the number of neutrons from 
induced and sponaneous fissions, respectively [15]. We note that for to — > 
— oo, the expectations of neutrons in the two regions ((Xi)) and ((N2)) will 
reach steady state due to the stationary source term with intensity S\. The 
solutions to the system of differential equations are found by differentiation 
of equation (J5J) with respect to (X, Y, Z) and then letting (X = Y = Z = 1). 
These read as 

{Ni) = Ri _ (-x 2 + u!x 2f )s iri 



AtiAt2 — (— Al + Z/^Ai/)( — A2 + Vi^2f) 

A 2 ffiri(A 2 -vf\2f) 
UJ1U2 

(N 2 ) = X 2 = ^\ (9) 
U1U2 

(Zt) = eXifNtf, (10) 
where e = A^/Ai/ and we have used the additional definitions u\ and u 2 as 

cot = -(A 2 + A 2 - (^Ai/ + z/iA 2 /)) 



+ 



^((Ai - A 2 ) - (i/^A/i - ^A 2/ )) 2 + 4A T1 A T2 , (11) 

w 2 = ^(A 2 + A 2 - (u^Xif + ^i 2 A 2/ )) 

- ^((Ai - A 2 ) - {v\\ fl - vi\2 f )) 2 + 4A T1 A T2 , (12) 

It can be noted that the expectation of the detections increases linearly with 
time and the number of neutrons in the second domain is directly determined 
by the transfer from the first to the second region A^i ■ The linear increase in 
the expectation of the detections is due to the fact that it is determined by 
integration of the expectation of the neutron number with respect to time. 
In order to find the variance of the detector counts we need to determine 
the second moments by yet another differentiation with respect to (X, Y, Z) 
followed by letting (X = Y = Z = 1). We find the variance of the detector 
counts through the relation cr| = (Zi) + fizz where the modified variance \izz 
is defined as \izz = (Z(Z — 1)) — (Z) 2 = o\ z — (Z) while in general we have 
Hxy = (XY) — (X)(Y). The differentiation procedure results in a system 
of six ordinary differential equations for the second order modified moments. 
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The differentiation procedure gives a system of six dynamical equations of 
the modified second moments as 

d 

— fixx = 2(-Ai + ulXif)fj, X x ~ 2\ T2 hxy + v 2 XifNt + Sir 2 , (13) 

-q^xy = (— Ai — A2 + v x Xif + v x X 2 f)/ixY — Xti^xx — AT 2 /iyy,(14) 
d 

-q^yy = 2(-A 2 + ^iA 2 /)/iyy - 2X T2 Hxy, (15) 

^/ixz = (-Ai + vl\if)tixz - X T2 ^yz + X d nxx, (16) 

= (-A 2 + z/i X 2 f)fj, Y z - Ari/ixz + A rf /i X y, (17) 

— /xzz = 2e\ lf nxz, (18) 

The equation system and its solution is rather analogous to the case of the 
Feynman-alpha equations in two-group theory and with one neutron energy 
group but including delayed neutrons as given in Ref. [13] and in Ref. [IB] , 
respectively Although an analytical solution for the general time-dependent 



system of equations (13) - (18) would be hard to find, we note that in the 
stationary state the system breaks down into two systems such that the 
solution of the first three equations is independent from the second, such that 
the moments (fJ>xx) = fixx, (Hxy) — fixY and (/iyy) = /iyy are constants. 
The equations describing detected particles need to be solved retaining the 
full time evolution by e.g. Laplace transforms. Moreover, it is found that the 
sought moment (fizz) is determined by quadrature of moment (fixz)- We 
find the constant 2nd moments as, 

((A 2 - ^i 2 A 2 /) 2 + wiu;2)(AT2A r 2 + 5ir 2 ) 

Vxx — — tj\ — n 2 r — \ > l iy J 

2(Ai + A 2 - u{Xif - u(\ 2f )ujiuj 2 

(A 2 - v\ X 2 f)X T \{X T2 N 2 + Sxr 2 ) 

2(Ai + A 2 - v{ Xif - v{X 2 f)u)iUJ 2 

X 2 T1 (X T2 N 2 + S ir2 ) 

^ YY 2(Ai + A 2 — vlXif — vfX 2 f)uiou 2 ' 

Here we have used the notation r 2 = d 2 r/dX 2 \x=i for the second factorial 



moment |15j . The objective now is to solve (16) and (17) by Laplace trans- 
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form methods and we find the transformed identity as, 

lSl\d^T2p>XY (s + X 2 — v\ X 2 f)\dfixX / 00 n 

^ XZ ~ sH(s) + JW(s) (22) 

with 

H(s) = s 2 + (co 2 + u l )s + u l u 2 . (23) 

Note that we have assumed that the initial values of the moments \xxz and 
[iyz were equal to zero at t — (at the start of the measurement), hence 
the roots of H(s) determine the temporal behavior of the Feynman-alpha 
formula. Moreover, the solution has many similarities to that found in Refs. 
[To*! [T3"] . We determine the variance to mean or Feynman-alpha formula 
by utilizing the relation for the variance a zz = (Z\) + [izz and after some 
algebra we find, 

a 2 (T) 1 - p-^ T 1 - P -^ T 

w- 1 ^ 1 "^^ 1 - 1 ^- (24) 

Here, the complete expressions for Y\ and Y 2 are quite lengthy, and are 
therefore given in the Appendix. However, it turns out that the sum Y = 
Yi + Y 2 takes a rather simple form that also determines the value of the 
Feynman-alpha for large gate times T — > 00 as, 

Y = Y 1 + Y 2 = u 2 XdiX2 -% X f XTlXT \ (25) 

In the next section we will consider some quantitative examples of the vari- 
ance to mean formula in Equation (24). 



3 Quantitative assessment of the variance to 
mean formula 

We will now discuss the quantitative application of the previously found 



variance to mean (Eq. 24) for a MOX fuel assembly, in DDAA and DDSI 
setups with varying amount of polyethylene moderator. In the simulations, 
the 58 fuel pins were recast into spherical geometry with the same amount 
of material, yielding a radius of roughly 13 cm. The spent fuel is of MOX 
type (U(94.14%)0 2 + Pu(5.86%)0 2 ) with approximately 4% fissile content. 
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We have considered two different thicknesses of the outer reflector shell with 
the sizes 10 and 500 cm where the larger moderator (500 cm) was assumed 
to be the equivalent of an infinite moderator. In the Monte-Carlo (MCNPX) 
simulations [I7J we consider two models, shown in Figure (II]). The first model 
represents a DDAA case (left), where the external pulsed neutron source 
was placed at the distance of 5 cm far from fuel sample having the same 
energy of neutrons as those produced in DD-fusion reactions. Note, that 
there is no speciel significance to the selection of DD-fusion reactions, the 
DT-fusion would have yielded similar results with the appropriate changes to 
the moderator thickness. The second model was constructed to be similar to 
a DDSI case (right), where the source was assumed to be equally distributed 
in the fuel volume and the energies of the source neutrons were sampled from 
the Watt fission spectrum, originating from the spontaneous fission of 240 Pu. 
The energy spectrum is determined by the function, 

p(E) ~ e ~ E/a smh(bE) 1/2 (26) 

where a = 0.799 MeV and b = 4.903 MeV^ 1 . Considering all possible cases 
we have hence carried out four MCNPX simulations in order to estimate the 
reaction intensities in two geometrical domains, fuel and moderator, with 
varying amount of moderation. In the simulations we have assumed a one 
energy group approximation. Reaction rates for each zone (and total) were 
determined by review of the neutron weight balance table (print table 130) 
of the MCNPX output which provides the fraction of source particles which 
enter each cell of the model, as it was done in papers of [H] and [19]. The 
simulation results are shown in Table 1 for the cases described above. 

Table 1: Simulation Results 



Intensity 

Ai [l/s] 
A 2 [l/s] 
Ati [l/s] 

AT2 [l/s] 

A/ [1/8] 



DDAA (500cm) 

0.7648 
1.5641 
0.5641 
0.4535 
0.1084 



DDAA (10cm) 

0.7441 
1.5502 
0.5501 
1.0185 
0.1049 



DDSI (500cm) 

1.6715 
2.0527 
1.2422 
0.8105 
0.2335 



DDSI (10cm) 

1.6665 
2.0453 
1.2406 
1.2588 
0.2317 



The simulation results are shown in Table[T]for the cases described above, 
where the reaction intensities obtained from the simulations were normalized 
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Figure 1: Show the geometry used for the Monte Carlo (MCNPX) simulations 
for the DDAA (left) and DDSI (right) cases in approximation of real and 
infinite moderator. 



Table 2: Results 



DDAA (500cm) DDAA (10cm) DDSI (500cm) DDSI (10cm) 
u eff 0.162 0.159 0.176 0.175 
wi [1/s] 0.6318 0.5299 0.9578 0.8360 
co 2 [1/s] 1.3837 1.4614 2.0916 2.1987 



to one starting neutron having a relative error less than 4%. In the in- 
terpretation of the simulation data, the sources in each region have to be 
properly accounted for remembering the specific setup of our analytical for- 
mula to achieve relevant comparisons. Withouth loss of generality in this 
study, we assume that there is no fission in the second region that contains 
only moderating material, by letting (Ai/ = Xf, v\ — v\ y A 2 / = = uf) and 
thus significantly simplifying the analytical formulas. In general we find a 
total balance of the intensities reflecting all the processes in the simulation 
(neglecting (n,xn) reactions that are several orders of magnitude smaller), 

1 + N fission = N ltf + N al + N a2 + N lost . (27) 

Here, on the left hand side the source is represented by 1 and Nf ission , where 
Nf ission is the number of fissions whereas the losses are represented by par- 
ticles lost in the fission process N^f signifying the number of fission events 
(the actual number of particles lost to fission) and the absorbtions in the 
two regions N a i and N a2 , respectively. Furthermore, since the moderating 
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material is finite, particles that may leave the system are accounted for by 
the term Ni ost . Note that all Ns are normalized to the number of source 
neutrons. Now we will consider the DDAA and DDSI situations separately 
as well as the two regions. In the DDSI situation we have a source in region 
I. For that region we have the balance, 

1 + N flssion = Ni t f + N al + N T r, (28) 

where Nt> = Nj^jj — Njj^j is the total flux of neutrons from region I to II 
and in region II we find 

N MI - Nj^j = N a2 + N lost . (29) 

In the DDAA case we find slightly different balance equations of the form in 
region I, 

N f ission - N T , = Ni t f + N al , (30) 
and in region II 

l + N T , = N a2 + N lost . (31) 

Here the source in region II is 1 + Nt>- In both the DDSI and DDAA cases, 
we identify Nt' with the net flux out of region I. Using the balance equations 
(27) - (31), we find the reaction intensities presented in Table [T] In order 



to check the consistency of our numerical studies we computed the effective 
multiplication factor v e ff- in the simulated cases this value can be found by 
using the definition, the ratio of number of neutrons to the neutrons in the 
preceding generation, 

»•» = m (32) 

where the probability of fission is the fission intensity divided by the total 
intensities \f/{\\ + A2). In combination with the multiplication factor u e ff 
we compute the time constant (wi,w 2 ) describing the system, and the values 
are given in Table [2] 

In order to evaluate analytical and numerical results for DDSI and DDAA 
cases, time distributions of the current of neutrons with energy higher than 
0.3 MeV crossing surface 1 (region I) in the direction of region II (tally Fl 
- MCNPX) were calculated. In Figure 2, one can see a comparison of the 



analytically estimated co± using Eq. (11) in the infinitely thick moderator 
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Figure 2: Fits for analytical depression for u x = 0.6318 tolhe DDAA "(left) 
and ui = 0.9578 to the DDSI (right). 




DDAA 500cm 
DDSI 500cm 



n ! i i i I 1 1 1 1 1 1 1 

o 50 ™ 150 200 10 20 30 40 50 

Figure 3: Evaluation of the variance to mean using the Values of the reaction 
intensities in Table. [I] for (left figure) DDAA, 500cm (blue line) and 10cm 
(red line) and DDSI 500cm (black line) and 10cm (magenta line); and for 
(right figure) DDAA, 500cm (blue line) and DDSI 500cm (black line). 



case and the simulation data of the neutron rates for DDAA (left) with 
lo\ = 0.6318 and DDSI (right) with oj\ = 0.9578. Good agreement is found 
between the numerical simulation and the analytical result. 

The variance to mean expression in Eq. ( 24 ) is displayed in Figure 3 using 
the data from Table 1 whereas the remaining parameters are ui = 2.80, 
z/ 2 = 4.635, \d = 0.1, ri = 1.0 and r 2 = 0.0. The results are shown for 
DDAA, 500cm (blue line) and 10cm (red line) and DDSI 500cm (black line) 
and 10cm (magenta line). The reaction intensities of the moderated cases are 
approximately of the same order of magnitude and thus the two exponentials 
are not easily distinguished. We find that there is a difference between the 
variance to mean for the DDAA and DDSI cases where the DDAA gives 
elevated values. For both DDAA and DDSI the infinitely moderated cases 
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have lower values compared to the cases with finite moderation. 



4 Discussion and conclusions 

We have developed a forward Kolmogorov approach for the two point theory 
of the Differential Die- A way Self Interrogation variance to mean (DSVM), 
including a compound Poisson source and the detection process. The results 
are an extension to those reported in [TJ El H3] where only one finite do- 
main was considered. However, we have restricted this study to be energy- 
independent and thus obtaining simplified analytical results. Retaining a 
two-point two-region theory, we expect to find up to four distinct exponen- 
tials. Furthermore, we have used Monte Carlo (MCNPX) simulations to find 
the reaction intensities needed to quantitatively assess the DSVM formula. 
We find that there is a difference between the variance to mean for the DDAA 
and DDSI cases where the DDAA gives elevated values in the limit of long 
gate times. For both DDAA and DDSI the infinitely moderated cases have 
lower values compared to the cases with finite moderation. In agreement 
with the results in [13j, we find that the presence of two exponents in the 
variance to mean (due to reaction intensities that are of the same order of 
magnitude) is most often not clearly visible. On one hand, this indicates that 
detection of the presence of fissile material may not be as obvious as with the 
Rossi-alpha method. On the other hand, two indistinguishable exponents in 
the variance to mean is not a problem in itself since determining uj\ and U2 
by curvefitting could be more accurate. Elucidating on the diagnostic value 
of the exponents in terms of determination of the sample parameters is not 
clear yet, and it requires further investigations, which will be reported in 
future work. 
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A The Fey nman- alpha formula 



In this appendix we give the full expressions of the Y\ and Y" 2 in Eq. (24), 

Yi = "2 j-i : yi, (1) 

v Arf(A 2 — ^1 A 2 /)AtiAt 2 w 2 , s 

y 2 — ^2 2~~2 1 2/2- (/J 

W 2 Wi + W 2 

Where yi and y 2 are found by, 

1 2(A 2 - z^ 2 A 2 /) 3 (u; 2 - Ai)(wi + w 2 ) - (AtiAt2) 2 w 2 



2/1 



+ 



2 (A 2 - z/ 2 A 2 /) 2 (a; 2 - wi)u;iu; 2 

AtiA T 2(A 2 - ^i 2 A 2/ )((Ai - v\\ lf )(u 2 - ^1) 
(A 2 - ^ 2 A 2 /) 2 (w 2 - Wi)o;ia; 2 



and 



V2 



(^1 + ^2)((A 2 - ^1 A 2/ ) + cji + 3u 2 )) 
(A 2 - ^ 2 A 2/ ) 2 (w 2 - u^u^t^ 



1 2(A 2 - ^iA 2/ ) 3 (u; 2 - (A 2 - z/ 2 A 2/ ))(u;i + w 2 ) + (A T iA T 2) u\ 

2 (A 2 - z/ 2 A 2 /) 2 (a; 2 - wi)u;iu; 2 
AnAy 2 (A 2 - z/ 2 A 2/ )((Ai - v\\xj){u2 - loi) 

(A 2 - ^1 A 2 /) 2 (w 2 - wi)wiw 2 
(wi + w 2 )((A 2 - v\\ 2 j) + 3a>i + cu 2 )) 
(A 2 - ^ 2 A 2/ ) 2 (cj 2 - w^wiu^ 



+ 



(4) 
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